######################################################
### OECD 2020 Risk That Matter ####
######################################################
# set the environment
library(haven);library(tidyverse);library(texreg)
library(xtable);library(psych);library(interflex)
# load the data 
rtm <- read_dta("rtm.dta")
#----------------------------------------------------------#
# dependent variable #
#----------------------------------------------------------#
# investment policy
rtm$invedu <- rtm$q32a
rtm$conben <- rtm$q32c
rtm$investment <- (rtm$q32a + rtm$q32b)
rtm$investment <- scale(rtm$investment)
# consumption policy 
rtm$consumption <- (rtm$q32c + rtm$q32d)
rtm$consumption <- scale(rtm$consumption)
#----------------------------------------------------------#
# independent variable #
#----------------------------------------------------------#
## automation threats ## 
# 27. how likely in five years? 
# a. My job will be replaced by a robot, computer software, an algorithm, or artificial intelligence.
rtm$threat1 <- as.numeric(rtm$q27a)
rtm$threat1[rtm$threat1 > 4] <- NA
# b. My job will be replaced by a person providing a similar service on an internet platform.
rtm$threat2 <- as.numeric(rtm$q27b)
rtm$threat2[rtm$threat2 > 4] <- NA
# c. I will lose my job because I am not good enough with new technology or because I will be replaced by someone with better technological skills.
rtm$threat3 <- as.numeric(rtm$q27c)
rtm$threat3[rtm$threat3 > 4] <- NA
# all: a+b+c 
rtm$threat <- as.numeric((rtm$threat1 + rtm$threat2 + rtm$threat3)/3)
# scale
rtm$threat1 <- scale(rtm$threat1)
rtm$threat2 <- scale(rtm$threat2)
rtm$threat3 <- scale(rtm$threat3)
rtm$threat <- scale(rtm$threat)
#----------------------------------------------------------#
# control variables #
#----------------------------------------------------------#
# s2 Sex
rtm$s2[rtm$s2 == 3] <- NA
rtm$female <- ifelse(rtm$s2 == 2, 1, 0)
# s3_age Age
rtm$age <- scale(rtm$s3_age)
# s6 Education
rtm$education <- scale(rtm$s6)
# s8_dec House income decile 1-10
rtm$income <- scale(rtm$s8_dec)
# college
rtm$college <- ifelse(rtm$s6 > 7, 1, 0)
# final data
rtm <- rtm %>%
  dplyr::select(investment, consumption, threat, threat1, threat2, threat3,
         education, income, female, age, college, ctrcode, invedu, conben)
write.csv(rtm, file = "rtm.csv")
#----------------------------------------------------------#
# main analysis #
#----------------------------------------------------------#
rtm <- read.csv("rtm.csv")
###############################
## Figure 1 ##
###############################
# X: IV country mean Y: DV country mean 
fig1 <- rtm %>%
  group_by(ctrcode) %>%
  mutate(Consumption = mean(consumption, na.rm = TRUE),
         Investment = mean(investment, na.rm = TRUE),
         Threat = mean(threat, na.rm = TRUE)) %>%
  distinct(Consumption, Investment, Threat)

fig1a <- ggplot(fig1, mapping = aes(x = Threat, y = Consumption)) +
  geom_text(aes(label = ctrcode)) +
  geom_smooth(method = "lm", col = "red") +
  xlab("Job Insecurity") +
  ylim(-0.3, 0.35) +
  xlim(-0.5, 0.8) +
  ylab("Preferences for Consumption Policy") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.y = element_text(size = 18))
fig1a
ggsave("figure1a.pdf")

fig1b <- ggplot(fig1, mapping = aes(x = Threat, y = Investment)) +
  geom_text(aes(label = ctrcode)) +
  geom_smooth(method = "lm", col = "red") +
  xlab("Job Insecurity") +
  ylab("Preferences for Investment Policy") +
  ylim(-0.4, 0.25) +
  xlim(-0.5, 0.8) +
  theme_bw() +
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.y = element_text(size = 18))
fig1b
ggsave("figure1b.pdf")

###############################
## Table 1 ##
###############################
# threat: all
inv1 <- lm(investment ~ threat + factor(ctrcode), data = rtm)
inv2 <- lm(investment ~ threat + education + income +
             female + age + factor(ctrcode), data = rtm)
con1 <- lm(consumption ~ threat + factor(ctrcode), data = rtm)
con2 <- lm(consumption ~ threat + education + income +
             female + age + factor(ctrcode), data = rtm)
texreg(list(inv1, inv2, con1, con2), 
          style = "AJPS", omit = "factor")

###############################
## Figure 2 ##
###############################
# consumption #  
# get Threat coefficients and standard error for coefficient plot
config <- lm(consumption ~ threat3 + education + income +
             female + age + factor(ctrcode), data = rtm)
summary(config) 

# data for figure 2a
fig2con <- data.frame(Threat = c("All", "Robot", "Internet", "Technology"),
                      pe = c(0.057488, 0.0538074,  0.0505419, 0.037564),
                      se = c(0.007086, 0.0067749, 0.0068001, 0.006764))
fig2con$lower <- fig2con$pe - (1.96*fig2con$se)
fig2con$upper <- fig2con$pe + (1.96*fig2con$se)

# figure 2a
fig2a <- ggplot(data=fig2con, 
                mapping=aes(x = Threat, y = pe, ymin = lower, ymax = upper, shape = Threat)) + 
  geom_pointrange(show.legend = FALSE, fatten = 9) + 
  ylim(-0.13, 0.13) +
  labs(y= 'Preferences for Social Consumption', x = 'Job insecurity: "My job will be replaced by..."') +
  theme_bw() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.y = element_text(size = 18))
fig2a
ggsave("figure2a.pdf")

## Investment ## 
invfig <- lm(investment ~ threat2 + education + income +
             female + age + factor(ctrcode), data = rtm)
summary(invfig)

# data for figure 2b
fig2inv <- data.frame(Threat = c("All", "Robot", "Internet", "Technology"),
                      pe = c(-0.092422, -0.076008, -0.074885, -0.070350),
                      se = c(0.007067, 0.006755, 0.006783, 0.006747))
fig2inv$lower <- fig2inv$pe - (1.96*fig2inv$se)
fig2inv$upper <- fig2inv$pe + (1.96*fig2inv$se)

fig2b <- ggplot(data = fig2inv, mapping=aes(x = Threat, y = pe, 
                                          ymin = lower, ymax = upper, shape = Threat)) + 
  geom_pointrange(show.legend = FALSE, fatten = 9) + 
  ylim(-0.13, 0.13) +
  labs(y= "Preferences for Social Investment", x = 'Job insecurity: "My job will be replaced by..."') +
  theme_bw() +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.y = element_text(size = 18))
fig2b
ggsave("figure2b.pdf")
###############################
## Appendix A.1&2 ##
###############################
### summary statistics ###
stat <- subset(rtm, select = c(investment, consumption, threat, threat1, threat2, threat3, female, education, income, age))
stat <- describe(stat)
stat <- subset(stat, select = c(n, mean, median, sd, min, max))
xtable(stat)
###############################
## Appendix A.3. ##
###############################
# threat1
inv1 <- lm(investment ~ threat1 + factor(ctrcode), data = rtm)
inv2 <- lm(investment ~ threat1 + education + income +
             female + age + factor(ctrcode), data = rtm)
con1 <- lm(consumption ~ threat1 + factor(ctrcode), data = rtm)
con2 <- lm(consumption ~ threat1  + education + income +
             female + age + factor(ctrcode), data = rtm)
texreg(list(inv1, inv2, con1, con2), style = "AJPS", 
          omit = "factor")
###############################
## Appendix A.4. ##
###############################
inv1 <- lm(investment ~ threat2 + factor(ctrcode), data = rtm)
inv2 <- lm(investment ~ threat2  + education + income +
             female + age + factor(ctrcode), data = rtm)
con1 <- lm(consumption ~ threat2 + factor(ctrcode), data = rtm)
con2 <- lm(consumption ~ threat2 + education + income +
             female + age + factor(ctrcode), data = rtm)
texreg(list(inv1, inv2, con1, con2), style = "AJPS", 
          omit = "factor")
###############################
## Appendix A.5. ##
###############################
inv1 <- lm(investment ~ threat3 + factor(ctrcode), data = rtm)
inv2 <- lm(investment ~ threat3  + education + income +
             female + age + factor(ctrcode), data = rtm)
con1 <- lm(consumption ~ threat3 + factor(ctrcode), data = rtm)
con2 <- lm(consumption ~ threat3 + education + income +
             female + age + factor(ctrcode), data = rtm)
texreg(list(inv1, inv2, con1, con2), style = "AJPS", 
          omit = "factor")
###############################
## Appendix A.6. ##
###############################
# skill level: college vs. non-college

c1 <- lm(investment ~ threat + college + factor(ctrcode), data = rtm)
c2 <- lm(investment ~ threat  + college + income +
             female + age + factor(ctrcode), data = rtm)
c3 <- lm(investment ~ threat*college + income +
           female + age + factor(ctrcode), data = rtm)
c4 <- lm(consumption ~ threat + college + factor(ctrcode), data = rtm)
c5 <- lm(consumption ~ threat + college + income +
             female + age + factor(ctrcode), data = rtm)
c6 <- lm(consumption ~ threat*college + income +
           female + age + factor(ctrcode), data = rtm)
texreg(list(c1, c2, c3, c4, c5, c6), style = "AJPS", 
          omit = "factor")

###############################
## Appendix Table 15 ##
###############################
rtm$invedu <- scale(rtm$invedu)
rtm$conben <- scale(rtm$conben)
inv1 <- lm(invedu ~ threat + factor(ctrcode), data = rtm)
inv2 <- lm(invedu ~ threat + education + income +
             female + age + factor(ctrcode), data = rtm)
con1 <- lm(conben ~ threat + factor(ctrcode), data = rtm)
con2 <- lm(conben ~ threat + education + income +
             female + age + factor(ctrcode), data = rtm)
texreg(list(inv1, inv2, con1, con2), 
       style = "AJPS", omit = "factor")

